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Abstract 

We  consider  a  variational  formulation  based  on  Maxwell’s  equations  for  the  prop¬ 
agation  of  high  frequency  (gigahertz  to  terahertz)  ultrashort  input  pulses  in  dielectric 
materials  modeled  by  a  linear  Debye  medium.  We  demonstrate  computationally  the 
emergence  of  Brillouin  precursors  in  the  material  (water)  and  the  fact  that  the  peak 
of  this  transient  is  attenuated  at  a  much  slower  rate  than  is  the  carrier  frequency.  In 
the  0.1  to  1  THz  regime  the  carrier  frequency  does  not  propagate  in  our  calculations. 
Only  the  precursors  enter  the  material,  and  this  is  in  line  with  experiments  reported 
by  Pleshko  and  Palocz  [11].  We  also  implement  models  that  include  nonlinearly  forced 
Debye  and  nonlinear  Debye  polarization  dynamics  and  demonstrate  the  importance 
of  nonlinear  effects,  especially  when  the  amplitude  of  the  input  signal  is  large.  This 
is  an  important  step  in  understanding  high  frequency  pulse  propagation,  and  it  has 
potential  applications  in  the  assessment  of  safety  standards  and  in  extending  current 
imaging  capabilities  in  both  civilian  and  military  uses. 


1  Introduction 

In  recent  years  there  has  been  a  great  technological  development  in  methodology  for 
generating  ultrashort  pulses  of  light  in  the  terahertz  range.  A  single  pulse  of  a  few 
picoseconds  long  allows  the  full  terahertz  spectrum  to  be  measured,  and  thus  facilitates 
better  (spatial  and  electromagnetic)  and  faster  imaging  results  than  more  traditional 
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techniques.  In  this  frequency  range  several  organic  materials  exhibit  strong  absorption 
and  dispersion  that  are  characteristic  to  the  particular  molecular  structure  and  content. 
Thus,  detection  of  the  temporal  distortions  of  the  electromagnetic  wave  caused  by  the 
interaction  yields  information  about  the  composition  of  the  material  in  real  time.  For 
example,  it  is  known  that  cancerous  and  benign  tumors  have  different  electromagnetic 
characteristics.  Thus  an  imaging  device  based  on  THz  waves  could  not  only  give 
information  about  the  structure  of  an  object  (geometrical  properties)  but  could  help 
in  determining  its  composition  and  electromagnetic  properties  as  well  in  a  non-invasive 
way. 

In  this  paper  we  present  extensive  numerical  simulations  modeling  the  propaga¬ 
tion  of  electromagnetic  waves  in  the  gigahertz  to  terahertz  frequency  range  through 
dielectric  materials  with  special  emphasis  on  ultrashort  input  pulses.  In  the  first  set 
of  numerical  experiments  we  concentrate  on  a  linear  Debye  medium,  and  study  the 
propagation  characteristics,  in  particular,  temporal  transients  such  as  the  Brillouin 
precursors  as  they  change  with  the  increasing  carrier  frequency  of  the  input  pulse.  Our 
results  are  in  agreement  with  theoretical  and  experimental  observations  [8,  11].  It  is 
known  that  most  materials  exhibit  nonlinear  characteristics  in  their  interaction  with 
the  electromagnetic  wave,  especially  when  the  input  signal  has  a  large  amplitude  [2]. 
Our  models  incorporate  different  nonlinear  polarization  mechanisms  which  influence 
the  propagation  characteristics.  In  particular,  we  study  a  nonlinearly  forced  linear 
Debye  model  [6],  as  well  as  a  mechanistically  nonlinear  Debye  model.  Our  computa¬ 
tional  framework  is  based  on  a  variational  formulation  of  Maxwell’s  equations.  We 
believe  that  this  framework  provides  a  means  to  capture  important  dynamic  phenom¬ 
ena  associated  with  short  input  pulses  (which  has  been  demonstrated  for  microwave 
propagation),  and  is  amenable  to  both  theoretical  and  computational  investigations. 
The  simulations  reported  here  form  an  important  step  towards  using  electromagnetic 
pulses  for  interrogation  of  unknown  materials  with  potential  applications  in  indus¬ 
trial  inspection,  security  screening,  medical  diagnostics  and  several  other  civilian  and 
military  problems. 

The  paper  is  organized  as  follows:  in  Section  2  we  present  the  general  variational 
framework  describing  the  physical  problem,  we  outline  the  numerical  methods  used  and 
present  the  results  for  a  linear  Debye  medium.  In  Section  3  we  include  nonlinearities 
in  the  form  of  a  nonlinearly  forced  Debye  polarization  mechanism,  while  in  Section  4 
we  present  our  initial  results  for  a  nonlinear  Debye  model. 


2  Variational  framework  and  linear  results 

In  this  section  we  describe  a  model  that  was  initially  developed  for  the  propagation 
of  microwaves  in  dielectric  materials  in  the  monograph  [3].  The  authors  studied  both 
the  direct  and  the  inverse  problem,  and  investigated  the  feasibility  of  using  windowed 
pulses  from  antenna  sources  to  interrogate  the  material.  They  developed  a  theoretical 
and  computational  framework  for  using  reflected  signals  to  identify  both  geometric  and 
electromagnetic  properties  of  an  object.  The  inverse  problems  were  solved  either  by 
use  of  a  supraconductive  backing  behind  the  object  or  by  use  of  acoustic  waves  from 
which  the  signal  is  reflected. 
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Our  focus  here  is  on  understanding  the  propagation  characteristics  of  pulses  with 
different  carrier  frequencies,  and  for  this  purpose  we  set  up  a  simplified  model.  However, 
we  note  that  this  simplified  formulation  can  be  readily  extended  to  treat  interfaces  and 
the  more  general  inverse  problems.  We  consider  as  in  [3]  an  infinite  (in  the  x  and 
y  direction)  slab  of  homogeneous  material  of  width  £  with  faces  parallel  to  the  xy 
plane.  However,  here  we  place  an  antenna  inside  the  material,  usually  in  the  middle. 
The  input  signal  is  a  planar  electromagnetic  wave  polarized  with  oscillations  in  the  xz 
plane  only.  The  electromagnetic  field  E  is  reduced  to  one  nontrivial  component,  in 
the  x  direction,  at  all  points  of  the  slab,  and  it  is  homogeneous  in  intensity  in  the  x 
and  y  directions.  The  electromagnetic  flux  density  D  and  the  polarization  P  inherit 
this  uniform  directional  property  from  E.  With  these  assumptions  Maxwell’s  equations 
yield 

/io£o pE  +  Ato P  +  HqcfE  —  E"  =  —no  js  t  >  0,  0  <  z  <  £,  (2.1) 

where  E  =  E(t,z )  and  P  =  P(t,z)  denote  the  x  components  of  the  electric  field  and 
the  polarization,  respectively,  and  D  =  £q£tE  +  P.  The  parameter  no  is  the  magnetic 
permeability,  £q  is  the  electric  permittivity  of  free  space  and  £r  is  the  relative  permittiv¬ 
ity  (see  [3]  for  discussions).  Here  Js  is  the  source  current  density,  while  the  conduction 
current  density,  Jc,  is  assumed  to  be  given  by  Ohm’s  law  as  Jc  =  aE.  It  is  argued 
in  [1,  3]  that  this  simple  form  of  Ohm’s  law  together  with  the  introduction  of  nonlo¬ 
cality  in  time  through  an  electric  susceptibility  kernel  (polarization  dynamics)  capture 
the  frequency  dependence  of  the  conductivity  of  dispersive  media.  In  this  section  we 
consider  a  linear  polarization  mechanism  proposed  by  Debye  [7]  to  model  the  behavior 
of  materials  with  a  permanent  dipole  moment,  e.g.,  water  (orientational  polarization). 
This  is  given  by 

tP  +  P  =  £0(£s  ~  £oo)E,  t>  0,  (2.2) 

where  £s  is  the  static  relative  permittivity,  while  =  er.  We  note  that  this  model 
can  be  realized  with  electric  susceptibility  kernel  g(t)  =  e~tp£o(£s  —  £oo)/t  in  the  con¬ 
volution  expression  P  =  g  *  E.  For  our  numerical  computations  we  use  approximately 
absorbing  boundary  conditions  at  the  ends  of  the  slab  to  prevent  large  reflections  of 
the  waves 

=  0,  (2.3) 

=  0,  (2.4) 

with  c2  =  •  The  approximate  nature  of  these  boundary  conditions  does  not  cause 

a  problem  for  us,  since  in  each  case  we  make  sure  that  the  computational  domain  is 
sufficiently  large,  so  that  we  do  not  see  any  reflections  in  our  time  frame.  We  complete 
the  system  with  initial  conditions 


c 


E-E' 
E  +  E' 


2=0 


z=t. 


E(0,z) 

=  0, 

0  <  z  <£, 

(2.5) 

E(0,z) 

=  0, 

0  <  z  <  £, 

(2.6) 

P(0,z) 

=  o, 

0  <  z  <£. 

(2.7) 
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We  define  V  =  iV1(0,f'),  H  =  L2(0,£)  and  write  (2.1)  and  (2.2)  in  weak  form  as 


{Hq£q erE,  <p)  +  (hq crE,  ip)  +  (P,  p>)  +  {Er,  ip')  +  ^-E(0)(p(0)  +  JplE(£)tp(£) 

c  c 

=  ~tM)(js,<p)  (2-8) 

(tP  +  P,  ip)  =  {£o(es  -  £oo)E,ip)  (2.9) 


for  all  ip  €  V.  Well-posedness  results  given  in  Chapter  3  of  [3]  can  be  used  to  guarantee 
the  usual  existence,  uniqueness  and  continuous  dependence  of  solutions  for  this  system. 

For  computational  purposes  we  scale  the  time  variable  by  a  factor  of  c  and  polar¬ 
ization  P  by  a  factor  of  1/eo,  i.e.,  we  let  E  =  E(ct ),  P  =  1/eo P(ct).  We  express  P  from 
(2.2)  and  substitute  it  into  (2.8).  The  new  equations  in  the  scaled  variables  are  (where 
for  the  sake  of  simplicity  of  notation  we  drop  the  overtildes) 

(erE,  ip)  +  ((r)cr  +  £dX )E,  ip)  +  (A 2P,  ip)  -  (£dX2E,  ip)  +  (E',  tp')  +  ^£(0)^(0) 
+y/TsE(£)ip(£)  =  -r)(  js,<p)  (2-10) 

(P  +  A  P,ip)  =  (£dX  E,tp),  (2-11) 


where  we  introduced  the  notation  £d  =  es  —  £oo,  A  =  X-  and  rj  =  ^ c .  We  discretize  the 
problem  in  the  space  variable  using  a  first  order  Galerkin  finite  element  approximation. 
We  divide  the  interval  [0,  £]  into  N  equal  subintervals  at  the  points  z  =  jh,  j  =  0, . . .  TV, 
where  h  =  £/N,  and  construct  standard  piecewise  linear  splines  4>f(z)-  The  finite  di¬ 
mensional  approximating  subspaces  to  V  will  be  taken  to  be  VN  =  j  cJ)q  ,  4>\  , . . . , 

Now  we  approximate  E(t,  z )  and  P(t,  z)  in  this  space  as 


N 

E(t,  z)  «  EN(t,  z)  =  J2  ef  (i)0f  (^)> 

3=0 

N 

P(t,z)  «  PN(t,z )  =  ^2pf{t)(j)f(z). 
3=0 


(2.12) 

(2.13) 


We  choose  the  space  of  test  functions  to  be  VN  also,  and  thus  we  find  that  (2.10)-(2.11) 
lead  to 


£rMe  +  ((rja  +  £dX)Al  +  p£^BD)e  +  A 2 Alp  +  (K  —  A 2£dAI)e  —  —pj ,  (2.14) 
Alp  =  —A  Mp  +  £dX  Ale  (2.15) 

where  e  =  (e^,  ef , . . . ,  e$)T,  p  =  (Po  ,Pi  ,  ■  ■  ■  ,Pn)T ,  J  =  (Jo,  ■  ■  ■ ,  Jn)T  ,  with 
Ji  =  (js,<j>i),  <  =  0, ...  TV,  and 


Alij  =  4>j)  =  [  4>i(z)4’j(z)dz, 

Jo 

Kij  =  (0'(/>')  =  ^(z)^ j(z)dz. 
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Here  BD  denotes  a  matrix  with  BDqq  =  BD^n  =  1,  while  all  other  elements  of  the 
matrix  are  0.  We  now  employ  a  second  order  central  difference  approximation  in  the 
time  variable  to  solve  (2.14).  We  take 


e{tn)  =  en 


£n+l  _  2en  +  en- 1 

At2 


e(tn)  =  en 


en+ 1  -  en_1 
2A  t 


and 


e{tn)  =  en 


en+l  +en  +  en- 1 
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Now  the  Galerkin  approximation  (2.14)  reduces  to  a  linear  equation  that  can  be  solved 
for  en+1  given  en,en_1  and  pn.  We  have  initial  condition  e°  =  0  and  we  approximate 
e1  ~  —r]^-js(0,z).  To  obtain  pn  we  solve  (2.15)  using  a  Crank-Nicholson  method. 
This  approximation  method  is  overall  0(h 2)  when  At  =  0(h).  We  also  note  that  the 
method  is  unconditionally  stable. 

In  the  following  we  describe  our  simulations  for  the  problem  above.  Our  physical 
constants  are  £o  =  8.854  •  10-12,  c  =  2.9980  •  108  m/sec ,  po  =  1.2566  •  10-6  We  carried 
out  all  the  calculations  assuming  that  the  dielectric  material  is  water  with  parameters 
er  =  5.5,  £ s  =  80.1,  r  =  8.1  •  10~12sec,  cr  =  10 ~5Ohm~1.  We  conducted  a  series 
of  numerical  experiments  with  windowed  input  signal  (generated  by  an  antenna  at 
z  =  i/2) 

Js  =  Asm3(ujt)x[o,tf]{t)ti{z  ~  i/2), 

with  amplitude  A  =  1,  frequencies  uj  =1  GHz,  10  GHz,  100  GHz  and  1  THz,  and 
final  input  time  tf  =  10/cu,  i.e. ,  we  used  5  complete  periods  of  the  input  signal.  We 
represent  the  results  two  different  ways:  (a)  we  took  snapshots  of  the  electric  field  in 
the  material  at  several  points  in  time  and  (b)  we  recorded  the  electric  field  at  particular 
points  in  the  material  from  time  t  =  0  to  some  t  =  tend.  The  value  of  tend  is  chosen  in 
each  case  such  that  the  input  signal  reaches  the  sides  of  the  material,  but  reflections  do 
not  propagate  backwards  significantly.  Figures  1  and  2  depict  these  simulations  for  the 
input  signal  with  carrier  frequency  1  GHz,  while  Figures  3  and  4  present  the  results 
for  the  signal  with  input  frequency  10  GHz. 

We  can  clearly  see  the  emergence  of  the  Brillouin  precursors  in  both  cases.  It  is 
important  to  observe  that  while  the  signal  is  attenuated  as  it  propagates  through  the 
material,  the  rate  of  attenuation  of  the  transient  (precursor)  is  slower  than  that  of  the 
carrier  frequency  (which  decays  exponentially).  This  fact  has  important  implications 
for  the  determination  of  safety  standards  since  the  precursor  may  penetrate  much 
deeper  with  more  energy  than  does  the  carrier  frequency  part  of  the  signal.  The  signal 
attenuation  curve  is  shown  in  Figure  5  for  the  1  GHz  case.  We  remark  that  these 
results  are  in  agreement  with  those  reported  in  [1],  where  the  authors  find  the  electric 
field  in  terms  of  a  Fourier  series  after  developing  the  input  pulse  train  in  a  series.  The 
drawback  of  that  method  is  that  it  does  not  readily  extend  to  higher  frequencies  and 
to  the  treatment  of  the  inverse  or  imaging  problem.  In  particular,  it  is  very  difficult  to 
treat  the  interaction  of  transmitted  and  reflected  waves. 

We  next  report  on  the  results  with  input  frequencies  100  GHz  (Figures  6  and  7) 
and  1  THz  (Figures  8  and  9). 
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Prop,  in  material  at  t=9.99375  ns 


Prop,  in  material  at  t=30.0188  ns 


Prop,  in  material  at  t=20.0063  ns 


Figure  1:  The  electric  field  in  the  material  at  t=10,  20,  30  and  40  ns,  1  GHz  example. 


Signal  at  z=1 .5  (meters) 


Signal  at  z=2.1  (meters) 


Signal  at  z=2.25  (meters)  Signal  at  z=2.5  (meters) 


Figure  2:  The  electric  field  recorded  at  the  antenna  and  at  a  distance  of  0.6,  0.75 
and  1  m  from  the  antenna,  1  GHz  example. 
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Prop,  in  material  at  t=0.331875  ns 


Prop,  in  material  at  t=0. 999375  ns 


Prop,  in  material  at  t=0. 665625  ns 


Figure  3:  The  electric  field  in  the  material  at  t=0.33,  0.66,  0.99  and  2.33  ns,  10 
GHz  example. 


Prop,  in  material  —  Signal  at  z=0.15  (meters) 
30 

2  20 
a> 

7/5  10 

o 
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2 

S  -io 
1-20 
-30 

0  1  2  3  4  5 


t  (ns) 


Signal  at  z=0.16  (meters) 


Signal  at  z=0.225  (meters)  Signal  at  z=0.25  (meters) 


Figure  4:  The  electric  field  recorded  at  the  antenna  and  at  a  distance  of  1  cm, 
7.5  cm  and  10  cm  from  the  antenna,  10  GHz  example. 
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Figure  5:  Attenuation  of  the  peak  of  the  transient  (dashed  line)  and  the  carrier 
frequency  (solid  line). 


Prop,  in  material  at  t=0. 04975  ns 
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Prop,  in  material  at  t=0. 149875  ns 


Prop,  in  material  at  t=0. 09981 25  ns 


Figure  6:  The  electric  field  in  the  material  at  t=0.049,  0.099,  0.14  and  0.25  ns,  100 
GHz  example. 


We  can  observe  that  the  carrier  frequency  does  not  penetrate  the  material  in  these 
cases.  It  is  only  the  transient  that  propagates  inside  the  material.  This  result  is 
expected  based  on  theoretical  considerations  [8,  p.313.]  and  experimental  observations 


Signal  at  z=0.016  (meters) 


Prop,  in  material  —  Signal  at  z=0.015  (meters) 


Signal  at  z=0.0225  (meters) 


Signal  at  z=0.025  (meters) 


Figure  7:  The  electric  field  recorded  at  the  antenna  and  at  a  distance  of  1  mm,  7.5 
mm  and  10  mm  from  the  antenna,  100  GHz  example. 


Prop,  in  material  at  t=0. 124531  ns 


Prop,  in  material  at  t=0. 249688  ns 


Prop,  in  material  at  t=0.374844  ns 


Prop,  in  material  at  t=0.5  ns 


Figure  8:  The  electric  field  in  the  material  at  t=0. 12,  0.24,  0.37  and  0.5  ns,  1  THz 
example. 
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Prop,  in  material  —  Signal  at  z=0.00375  (meters)  Signal  at  z=0.005  (meters) 


Signal  at  z=0. 00625  (meters) 


Figure  9:  The  electric  field  recorded  at  the  antenna  and  at  a  distance  of  1.25  mm, 
2.5  mm  and  2.81  mm  from  the  antenna,  1  THz  example. 


[U]- 

Thus  we  can  see  that  our  numerical  simulations  based  on  a  variational  formulation 
provide  a  detailed  description  of  the  propagation  of  the  electromagnetic  pulse,  in  par¬ 
ticular,  we  are  readily  able  to  capture  the  temporal  transients,  such  as  the  Brillouin 
precursors. 


3  Nonlinear ly  forced  Debye  polarization  model 

In  this  section  we  introduce  a  polarization  dynamics  driven  by  a  nonlinear  function 
of  the  electric  field.  It  is  known  that  most  materials  exhibit  nonlinear  polarization 
characteristics  especially  when  the  amplitude  of  the  input  signal  is  large.  In  a  first 
attempt  to  incorporate  nonlinear  effects  in  the  Debye  polarization  model  we  considered 

tP  +  P  =  soEdE  +  f{E).  (3-1) 

We  showed  that  this  constitutive  relationship  together  with  the  weak  formulation  of 
the  problem  (2.8)  is  well-posed.  In  particular,  we  have  the  following  result 

Theorem  3.1  If  the  nonlinear  function  f  :  1R  — >  1R '.is  C1,  with  /( 0)  =  0,  and  f(z)  > 
0  for  all  z  £  IR,  and  f  is  affine  at  infinity,  i.e.,  there  exist  constants  R,L  such  that 
for  every  |x|  >  R ,  |/(x)|  <  L |x|,  then  there  exists  a  weak  solution  to  the  system  (2.8), 
(3.1)  with  initial  conditions 


E(0,z)  =  $, 


(3.2) 
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E(0,z)  =  *,  (3-3) 

for  any  $  E  V,  'F  £  H,  Js  E  H2(0,T\V*).  The  weak  solution  is  unique,  and  depends 
continuously  on  4>,  'P  and  Js  under  the  additional  assumption  that  f  is  globally  Lips- 
chitz. 

The  details  of  the  proofs  can  be  found  in  [6] . 

In  our  numerical  simulations  we  take  f(E)  =  /3_ZT3 ,  which  is  an  approximation  to 
a  saturation  limited  nonlinearity  required  by  the  above  theorem.  We  implemented 
the  same  numerical  method  outlined  in  the  previous  section,  except  that  in  this  case 
(2.14)-(2.15)  contain  nonlinear  terms.  We  used  a  functional  iteration  in  the  nonlinear 
version  of  (2.14)  to  obtain  en,  which  is  then  used  to  update  the  polarization  term 
pn.  Comparison  of  the  nonlinearly  forced  model  (3.1)  with  the  corresponding  linear 
dynamics  is  depicted  in  Figures  11  and  12  for  the  10  GHz  and  1  THz  case,  respectively. 
In  these  simulations  we  have  a  weak  nonlinearity  with  (3  =  —5  x  10-6  and  the 


Prop,  in  material  —  Signal  at  z=0.15  (meters) 


Signal  at  z=0.16  (meters) 


Signal  at  z=0. 15375  (meters) 


Figure  10:  Comparison  of  the  linear  (solid  line)  and  nonlinearly  forced  (dashed  line) 
results,  10  GHz  example. 


amplitude  of  the  input  signal  is  small,  A  =  10.  The  linear  and  nonlinear  results  do  not 
differ  substantially  in  the  10  GHz  case  (see  Figure  11).  In  that  case  if  (3  is  positive 
we  observed  that  a  nonlinearly  forced  polarization  dynamics  leads  to  a  signal  whose 
main  part  arrives  slightly  earlier  and  is  slightly  larger  than  the  corresponding  portion 
of  the  signal  in  the  linear  material.  The  reverse  is  true  if  the  E 3  term  has  a  negative 
coefficient  in  (3.1).  The  same  observation  holds  true  for  the  case  of  input  signals  with 
carrier  frequency  1  GHz.  However,  in  the  1  THz  example  (see  Figure  12),  the  situation 
is  different.  The  main  part  of  the  signal  arrives  at  approximately  the  same  time  in 
the  linear  and  nonlinearly  forced  cases,  and  we  can  see  a  larger  difference  between  the 
linear  and  nonlinearly  forced  results. 
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Signal  at  0.0025  (meters) 


Signal  at  z=0.00375  (meters) 


Signal  at  z=0.003  (meters) 


Figure  11:  Comparison  of  the  linear  (solid  line)  and  nonlinear ly  forced  (dashed  line) 
results,  1  THz  example. 


4  Nonlinear  Debye  Polarization  model 

Our  last  set  of  numerical  simulations  are  based  on  a  nonlinear  Debye  medium  with 
polarization  dynamics  given  by: 

rP  +  g(P)  =  £0(es  -  £oo)E.  (4.1) 

We  choose  g(P)  =  P  +  sP 3,  where  s  is  a  small  parameter.  The  cubic  nonlinear  term 
is  chosen  since  most  biological  tissues  are  not  expected  to  have  an  axis  of  symmetry. 
However,  it  must  be  considered  as  an  approximation  to  a  saturation  limited  nonlinear 
mechanism  that  can  be  anticipated  for  orientational  polarization.  Large  departures 
from  the  linear  behavior  can  be  observed  even  for  weak  nonlinearities  if  the  amplitude 
of  the  input  signal  is  sufficiently  large.  In  our  simulations  s  =  10~3  and  the  amplitude 
A  of  the  input  signal  is  in  the  range  1010  to  1012.  Figures  13  and  14  show  the  comparison 
between  the  nonlinear  Debye  model  (4.1)  and  the  corresponding  linear  dynamics  in  the 
10  GHz  and  1  THz  case,  respectively. 

We  make  several  interesting  observations.  With  A  =  1010  in  the  10  GHz  example 
(Figure  13),  it  is  clear  that  the  leading  edge  of  the  main  part  of  the  signal  (i.e. ,  the 
Brillouin  precursor)  arrives  faster  than  in  the  corresponding  linear  simulation.  More¬ 
over,  the  electric  field  is  considerably  larger  in  the  nonlinear  Debye  medium.  However, 
an  input  signal  with  amplitude  A  =  1010  leads  to  no  difference  in  the  electric  field 
in  the  linear  and  nonlinear  simulations  in  the  1  THz  example  (not  shown)  and  one 
must  increase  the  amplitude  to  see  any  difference  in  the  results  (see  Figure  14,  where 
A  =  1011).  The  difference  is  similar  to  that  we  observed  in  the  nonlinearly  forced  case 
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Figure  12:  Comparison  of  the  linear  (solid  line)  and  nonlinear  (dashed  line)  results, 
10  GHz  example. 


in  the  terahertz  regime.  The  main  part  of  the  signal  appears  at  the  same  time,  but  in 
the  nonlinear  material  the  field  is  much  larger,  almost  twice  as  large  as  in  the  linear 
case.  We  note  that  for  A  =  1012  the  field  produced  in  the  nonlinear  case  is  two  orders 
of  magnitude  larger  than  that  in  the  corresponding  linear  case. 


5  Conclusion 

In  this  paper  we  presented  numerical  results  for  the  propagation  of  high  frequency 
windowed  pulses  in  dielectric  materials.  Although  the  problem  was  simplified  by  gen¬ 
erating  the  signal  inside  the  material  we  emphasize  that  this  formulation  can  readily 
be  extended  to  treat  interfaces  and  reflected  signals,  and  is  therefore  suitable  to  treat 
the  inverse  or  interrogation  problem.  The  understanding  and  correct  description  of 
temporal  transients  is  important  in  its  own  right  in  assessing  safety  standards  and 
making  the  interrogation  problem  practically  feasible.  As  we  noted  in  the  Introduc¬ 
tion,  electromagnetic  interrogation  with  high  frequency  windowed  pulses  have  many 
potential  applications.  We  remark  that  it  is  desirable  to  extend  the  physical  problem  to 
two  and  three  dimensions,  to  create  more  realistic  models  for  these  applications.  The 
two  dimensional  extension  has  been  carried  out  for  a  linear  Debye  medium  in  [5]  using 
perfectly  matched  layers.  It  was  also  implemented  under  slightly  different  assumptions 
and  using  a  different  numerical  method  in  [4]. 
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x  1 09  Signal  at  z=0. 0035  m 


Figure  13:  Comparison  of  the  linear  (solid  line)  and  nonlinear  (dashed  line)  results, 
1  THz  example. 
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